function [TP_a,tp_a] = AES2019P5(rho,i,I,T,L,Draw,Progress,c,d,x,vc0,sigmasq,alphac,alphad,alphax,alphasigma,beta,gamma,mu,nux,nusigma,sigmabar,psic,psid,chi,gap_x,gap_sigmasq,Grid_x,Grid_sigmasq)
%AES2019P5 Calculate the timing premium in the horizon-dependent risk
%aversion model (Andreis, Eisenbach, and Schmalz 2019) when rho is
%different from 1
%   Input: rho controls the intertemporal rate of substitution; T, L, Draw,
%   and Progress are global settings imported from AES2019_2; all other
%   variables are the same as in AES2019_8. Output: TP_a is the timing
%   premium, and tp_a is the (negative) log timing premium.

% Initialisation
vc1 = zeros(T,1); % next period vc
v1_d = zeros(T,1); % value with early resolution

% Expectation for every period
for t=1:T

    % Step 1: with uncertainty
    % vc ratio
    f_x = min(floor((x(t)+20.5*gap_x)/gap_x),40); % position, max 40                 
    f_sigmasq = min(floor((sigmasq(t)-sigmabar^2+12*gap_sigmasq)/gap_sigmasq),40); % max 40                   
    if f_x>0 && f_sigmasq>0
        d_x = (x(t)-Grid_x(f_x,1))/gap_x; % distance from the grid (0 to 1)
        d_sigmasq = (sigmasq(t)-Grid_sigmasq(f_sigmasq))/gap_sigmasq;
        vc1(t) = vc0(f_x,f_sigmasq)+(vc0(f_x+1,f_sigmasq)-vc0(f_x,f_sigmasq))*(1-d_x)+(vc0(f_x+1,f_sigmasq+1)-vc0(f_x+1,f_sigmasq))*(1-d_sigmasq);   
    elseif f_x<=0 && f_sigmasq>0
        d_sigmasq = (sigmasq(t)-Grid_sigmasq(f_sigmasq))/gap_sigmasq;
        vc1(t) = vc0(1,f_sigmasq)+(vc0(1,f_sigmasq+1)-vc0(1,f_sigmasq))*(1-d_sigmasq); 
    elseif f_x>0 && f_sigmasq<=0
        d_x = (x(t)-Grid_x(f_x,1))/gap_x; % distance from the grid (0 to 1)
        vc1(t) = vc0(f_x,1)+(vc0(f_x+1,1)-vc0(f_x,1))*(1-d_x); 
    else
        vc1(t) = vc0(1,1);
    end

    % Step 2: early resolution
    v2_d = zeros(Draw,1); % realisation of second-period value

    % Initialisation
    for draw=1:Draw
    sigmasq0 = sigmasq(t); % conditional volatility of consumption
    x0 = x(t); % varying component of the expected consumption growth rate 
    d0 = d(t); % initial log dividend (without loss of generality)
    c0 = c(t); % initial log consumption (without loss of generality)
    k = (0:L-1)'; % 0 to L-1 vector

    % Empty vectors
    sigmasq_d = zeros(L,1); % conditional volatility of consumption
    x_d = zeros(L,1); % varying component of the expected consumption growth rate 
    d_d = zeros(L,1); % log dividend
    c_d = zeros(L,1); % log consumption

    % Random draws
    epsilon_d = randn(L,4);
 
    % Volatility
    sigmasq_d(1) = sigmabar^2+nusigma*(sigmasq0-sigmabar^2)+alphasigma*epsilon_d(1,1); % first value (period t+1)
    for l=2:L
        sigmasq_d(l) = sigmabar^2+nusigma*(sigmasq_d(l-1)-sigmabar^2)+alphasigma*epsilon_d(l,1);
        if sigmasq_d(l)<sigmabar^2/10000
            sigmasq_d(l) = sigmabar^2/10000;
        end
    end

    % Varying component of the expected consumption growth rate 
    x_d(1) = nux*x0+alphax*(sigmasq0^0.5)*epsilon_d(1,2); % first value (period t+1)
    for l=2:L
        x_d(l) = nux*x_d(l-1)+alphax*(sigmasq_d(l-1)^0.5)*epsilon_d(l,2);
    end

    % Dividend
    d_d(1) = d0+mu+psic*x0+alphad*(sigmasq0^0.5)*epsilon_d(1,3)+chi*alphac*(sigmasq0^0.5)*epsilon_d(1,4);
    for l=2:L
        d_d(l) = d_d(l-1)+mu+psid*x_d(l-1)+alphad*(sigmasq_d(l-1)^0.5)*epsilon_d(l,3)+chi*alphac*(sigmasq_d(l-1)^0.5)*epsilon_d(l,4);
    end

    % Consumption
    c_d(1) = c0+mu+psic*x0+alphac*(sigmasq0^0.5)*epsilon_d(1,4); % first value (period t+1)
    for l=2:L
        c_d(l) = c_d(l-1)+mu+psic*x_d(l-1)+alphac*(sigmasq_d(l-1)^0.5)*epsilon_d(l,4);
    end

    % Second period
    v2_d(draw) = (1/(1-rho))*log(((1-beta)*(beta.^k))'*exp((1-rho)*c_d)); % sum of future consumptions (in one draw)
    end

    % First period
    v1_d(t) = real((1/(1-rho))*log((1-beta)*exp((1-rho)*c(t))+beta*mean(exp(v2_d).^(1-gamma))^((1-rho)/(1-gamma))));

    % Progress
    if Progress==1
        if T*I<=1000
            fprintf('Progress: %4.1f%s\n',100*((i-1)*T+t)/(I*T),'%')
        else
            progressbar = ceil(1:(T-1)/(floor(1000/I)-1):T);
            if any(progressbar==t)
                fprintf('Progress: %4.1f%s\n',100*((i-1)*T+t)/(I*T),'%')
            end
        end
    end
end

% Step 3: time permium 
v1 = vc1+c; % value with uncertaiinty
tp_a = v1-v1_d; % alternative time premium
TP_a = mean(ones(T,1)-exp(tp_a));

% End of function
end